###################################################
######                                       ######
###### BEN RADFORD                           ######
###### POLSCI 303                            ######
###### MAXIMUM LIKELIHOOD ESTIMATION         ######
###### JOCR REPLICATION                      ######
###### Fuhrmann & Kreps                      ######
###### "Targeting Nuclear Programs..."       ######
######                                       ######
###################################################


#############################
#                           #
#  CHANGE DEPENDING ON PC:  #
#                           #
#############################
rm(list=ls(all=TRUE)) 
#setwd("C:\\Users\\Ben\\Documents\\School\\")
setwd("A:\\Spring 2012\\PS303 - MLE\\JOCR Rep\\")
#setwd("F:\\")


######################
#                    #
#  INSTALL PACKAGES  #
#                    #
###################### 

install.packages("xtable")			#INSTALL LATEX EXPORT PACKAGE
library(xtable)
install.packages("apsrtable")			#INSTALL APSRTABLE FOR LATEX EXPORT
library(apsrtable)
install.packages("foreign")			#IMPORT DATA
library(foreign)	
install.packages("ggplot2")			#COOL GRAPHICS
library(ggplot2)
install.packages("lattice")			#MORE COOL GRAPHICS
library(lattice)
install.packages("gmodels")			#FOR CROSSTABLE
library(gmodels)
install.packages("sandwich")			#RUBUST STANDARD ERRORS
library(sandwich)
install.packages("car")				#FOR HCCM
library(car)
install.packages("Zelig")			#FOR SOME COOL REGRESSION STUFF
library(Zelig)
install.packages("lmtest")
library(lmtest)
#install.packages("arm")			#DISABLED FOR CONFLICT WITH sim IN ZELIG
#library(arm)					#DISABLED FOR CONFLICT WITH sim IN ZELIG
install.packages("QuantPsyc")			#STANDARDIZED COEFFICIENTS IN LM
library(QuantPsyc)
install.packages("mass")
library(MASS)
install.packages("separationplot")
library('separationplot')


#############################
#                           #
# LOAD DATA FOR REPLICATION #
#                           #
#############################

#repdata <- read.dta("Targeting Nuclear Programs.dta")
#head(repdata)
#cowcodes <- read.csv("cowcodes.csv")

load("datawithnames.RData")

head(repdata)

length(repdata[,1])-sum(as.numeric(complete.cases(repdata)))

##################################
#                                #
# SAVE NEW DATA W/ COUNTRY NAMES #
#                                #
##################################
#repdata <- data.frame(repdata)

#for (i in 1:length(repdata$ccode1)){
#  repdata$cname1[i] <- as.character(cowcodes$StateAbb[repdata$ccode1[i]==cowcodes$CCode])
#  repdata$cname2[i] <- as.character(cowcodes$StateAbb[repdata$ccode2[i]==cowcodes$CCode])
#  print(i/length(repdata$ccode1))
#}

#save(repdata,file="datawithnames.RData")

###################################
# CREATE CROSS TABS/SUMMARY STATS #
###################################
# hostileMID  = Conflict Y/N?     #
# polity2DUM  = Polity score      #
# s_un_regDUM = For.Pol.Simi.     #
###################################
con.conf <- table(repdata[repdata$polrel==1,]$consider1,repdata[repdata$polrel==1,]$hostileMID)
att.conf <- table(repdata[repdata$polrel==1,]$attack1,repdata[repdata$polrel==1,]$hostileMID)
con.poli <- table(repdata[repdata$polrel==1,]$consider1,repdata[repdata$polrel==1,]$polity2DUM)
att.poli <- table(repdata[repdata$polrel==1,]$attack1,repdata[repdata$polrel==1,]$polity2DUM)
con.simi <- table(repdata[repdata$polrel==1,]$consider1,repdata[repdata$polrel==1,]$s_un_regDUM)
att.simi <- table(repdata[repdata$polrel==1,]$attack1,repdata[repdata$polrel==1,]$s_un_regDUM)

##############
# CHI2 TESTS #
##############
chisq.test(con.conf)
chisq.test(att.conf)
chisq.test(con.poli)
chisq.test(att.poli)
chisq.test(con.simi)
chisq.test(att.simi)

###############
# PROP TABLES #
###############
con.conf
att.conf
con.poli
att.poli
con.simi
att.simi
prop.table(con.conf,2)
prop.table(att.conf,2)
prop.table(con.poli,2)
prop.table(att.poli,2)
prop.table(con.simi,2)
prop.table(att.simi,2)

########################################
# HOW MANY ATTACKS/CONSIDERED ATTACKS? #
########################################
sum(repdata$attack1,na.rm=T)
sum(repdata$consider1,na.rm=T)

############################
# NORWAY 1941-1944 MISSING #
############################
repdata$year[repdata$cname1=="NOR" & repdata$cname2=="GMY"]
min(repdata$year[repdata$cname1=="NOR"])

repdata$year[repdata$cname1=="ISR" & repdata$cname2=="PAK"]
names(repdata)

###############################
#                             #
# TABLE 4 - RARE EVENTS LOGIT #
#                             #
###############################
# TABLE 4 MODEL 1
z.out1 <-zelig(data=repdata[repdata$polrel==1,], consider1 ~ hostileMID + polity2 + s_un_reg + noconsideryrs + X_spline1 + X_spline2 + X_spline3, model="relogit")
# TABLE 4 MODEL 2
z.out2 <-zelig(data=repdata[repdata$polrel==1,], consider1 ~ hostileMID + polity2 + s_un_reg + cap_2 + pwrratio + postArt56 + postCW + prgmyrs + NCAnosafetodate + cap_1 + contig + noconsideryrs + X_spline1 + X_spline2 + X_spline3, model="relogit")
# TABLE 4 MODEL 3
z.out3 <-zelig(data=repdata[repdata$polrel==1,], consider1 ~ hostileMID + postArt56 + postCW + prgmyrs + NCAnosafetodate + noconsideryrs + X_spline1 + X_spline2 + X_spline3, model="relogit")
# TABLE 4 MODEL 4
z.out4 <-zelig(data=repdata[repdata$polrel==1,], attack1 ~ hostileMID + polity2 + s_un_reg +  noattackyrs + X_spline1a + X_spline2a + X_spline3a, model="relogit")
# TABLE 4 MODEL 5
z.out5 <-zelig(data=repdata[repdata$polrel==1,], attack1 ~ hostileMID + polity2 + s_un_reg +  cap_2 + pwrratio + postArt56 + postCW + prgmyrs + NCAnosafetodate + cap_1 + contig + noattackyrs + X_spline1a + X_spline2a + X_spline3a, model="relogit")
# TABLE 4 MODEL 6
z.out6 <-zelig(data=repdata[repdata$polrel==1,], attack1 ~ hostileMID + postArt56 + postCW + prgmyrs + NCAnosafetodate + noattackyrs + X_spline1a + X_spline2a + X_spline3a, model="relogit")
summary(z.out1)
summary(z.out2)
summary(z.out3)
summary(z.out4)
summary(z.out5)
summary(z.out6)

apsrtable(z.out1,z.out2,z.out3,z.out4,z.out5,z.out6)

pred1 <- predict(z.out1,newdata = data.frame(repdata))
pred1 <- data.frame(pred1)
pred1 <- cbind(pred1, repdata$attack1)
pred1 <- pred1[complete.cases(pred1),]
separationplot(pred=pred1[,1], pred1[,2])

pred2 <- predict(z.out2,newdata = data.frame(repdata))
pred2 <- data.frame(pred2)
pred2 <- cbind(pred2, repdata$attack1)
pred2 <- pred2[complete.cases(pred2),]
separationplot(pred=pred2[,1], pred2[,2])

pred3 <- predict(z.out3,newdata = data.frame(repdata))
pred3 <- data.frame(pred3)
pred3 <- cbind(pred3, repdata$attack1)
pred3 <- pred3[complete.cases(pred3),]
separationplot(pred=pred3[,1], pred3[,2])

pred4 <- predict(z.out4,newdata = data.frame(repdata))
pred4 <- data.frame(pred4)
pred4 <- cbind(pred4, repdata$attack1)
pred4 <- pred4[complete.cases(pred4),]
separationplot(pred=pred4[,1], pred4[,2])

pred5 <- predict(z.out5,newdata = data.frame(repdata))
pred5 <- data.frame(pred5)
pred5 <- cbind(pred5, repdata$attack1)
pred5 <- pred5[complete.cases(pred5),]
separationplot(pred=pred5[,1], pred5[,2])

pred6 <- predict(z.out6,newdata = data.frame(repdata))
pred6 <- data.frame(pred6)
pred6 <- cbind(pred6, repdata$attack1)
pred6 <- pred6[complete.cases(pred6),]
separationplot(pred=pred6[,1], pred6[,2])

##########################################
# FIRST DIFFERENCES PLOT TABLE 4 MODEL 2 #
##########################################
minMids<-c(1, 0,summary(repdata$polity2)[1],0,summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noconsideryrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
maxMids<-c(1, 0,summary(repdata$polity2)[6],0,summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noconsideryrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
outcomeMax<-matrix(NA, nrow=1000,ncol=1)
outcomeMin<-matrix(NA, nrow=1000,ncol=1)

#draw1<-mvrnorm(1,z.out2$coef,vcov(z.out2))
for (i in 1:1000) 
{
  outcomeMax[i]<-as.vector(exp(draw1%*%maxMids)/(1+exp(draw1%*%maxMids)))
  outcomeMin[i]<-as.vector(exp(draw1%*%minMids)/(1+exp(draw1%*%minMids)))

}

plot(density((outcomeMax)),xlab="Probability of Considered Attack",ylab="",xlim=c(0,.02),ylim=c(0,300), las=1,main="Probability of Attack\n Given Target's Polity",bty="n",lwd=3,col="blue")
lines(density((outcomeMin)),lwd=3,col="red")
legend(.01,300,c("Polity 2 score -9","Polity 2 score 10"), lwd=c(3,3), lty=c(1,1), col=c("Blue","Red"))

##########################################
# FIRST DIFFERENCES PLOT TABLE 4 MODEL 5 #
##########################################
minMids<-c(1, 0,summary(repdata$polity2)[4],1,summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noattackyrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
maxMids<-c(1, 1,summary(repdata$polity2)[4],1,summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noattackyrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
outcomeMax<-matrix(NA, nrow=1000,ncol=1)
outcomeMin<-matrix(NA, nrow=1000,ncol=1)

#draw1<-mvrnorm(1,z.out5$coef,vcov(z.out5))
for (i in 1:1000) 
{
  outcomeMax[i]<-as.vector(exp(draw1%*%maxMids)/(1+exp(draw1%*%maxMids)))
  outcomeMin[i]<-as.vector(exp(draw1%*%minMids)/(1+exp(draw1%*%minMids)))
}

plot(density((outcomeMax)),xlab="Probability of Attack",ylab="",xlim=c(0,.015),ylim=c(0,2000), las=1,main="Probability of Attack\n Given Prior Militarized Conflict",bty="n",lwd=3,col="blue")
lines(density((outcomeMin)),lwd=3,col="red")
legend(.005,2000,c("Prior militarized conflict","No prior militarized conflict"), lwd=c(3,3), lty=c(1,1), col=c("Blue","Red"))



#################################
# BETTER FIRST DIFFERENCES PLOT #
#################################
draw1 <- mvrnorm(1000,z.out2$coefficients,vcov(z.out2))
sim1 <- c(1, 0, -10 ,summary(repdata$s_un_reg)[1],summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noattackyrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
final <- (exp(draw1%*%sim1))/(1+exp(draw1%*%sim1))
sim2 <- c(1, 1, -10, summary(repdata$s_un_reg)[1],summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noattackyrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
final2 <- (exp(draw1%*%sim2))/(1+exp(draw1%*%sim2))

for (i in -9:10){
sim1 <- c(1, 0,i,summary(repdata$s_un_reg)[1],summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noattackyrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
pis1 <- (exp(draw1%*%sim1))/(1+exp(draw1%*%sim1))
final <- cbind(final,pis1)
}
for (i in -9:10){
sim2 <- c(1, 1,i,summary(repdata$s_un_reg)[1],summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noattackyrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
pis2 <- (exp(draw1%*%sim2))/(1+exp(draw1%*%sim2))
final2 <- cbind(final2,pis2)
}

for (i in 1:21){
final[,i] <- sort(final[,i])
}
for (i in 1:21){
final2[,i] <- sort(final2[,i])
}

plot(c(),c(),xlim=c(-10,10),ylim=c(0,.12),main="Systematic Uncertainty in Predicted Probabilities\n of Considered Attack (95% CI)", ylab="Probability of Considered Attack", xlab="Target Polity Score")
for (i in 1:21){
lines(x=c(i-11,i-11),y=c(final[25,i],final[975,i]),col="Blue",lwd=2)
lines(x=c(i-10.9,i-10.9),y=c(final2[25,i],final2[975,i]),col="Red",lwd=2)
}
legend(0,.12,c("Prior militarized conflict","No prior militarized conflict"),lty=c(1,1),lwd=c(2,2),col=c("Red","Blue"))




##############################
#                            #
#   TABLE 4 ORDERED LOGIT    #
#                            #
##############################
repdata$DVordered <- (repdata$consider1-repdata$attack1)+(2*repdata$attack1)

z.order1 <- zelig(data=repdata[repdata$polrel==1,], as.ordered(DVordered) ~ hostileMID + polity2 + s_un_reg + noconsideryrs + X_spline1 + X_spline2 + X_spline3, model="ologit")
# TABLE 4 MODEL 2
z.order2 <- zelig(data=repdata[repdata$polrel==1,], as.ordered(DVordered) ~ hostileMID + polity2 + s_un_reg + cap_2 + pwrratio + postArt56 + postCW + prgmyrs + NCAnosafetodate + cap_1 + contig + noconsideryrs + X_spline1 + X_spline2 + X_spline3, model="ologit")
# TABLE 4 MODEL 3
z.order3 <- zelig(data=repdata[repdata$polrel==1,], as.ordered(DVordered) ~ hostileMID + postArt56 + postCW + prgmyrs + NCAnosafetodate + noconsideryrs + X_spline1 + X_spline2 + X_spline3, model="ologit")
# TABLE 4 MODEL 4
z.order4 <- zelig(data=repdata[repdata$polrel==1,], as.ordered(DVordered) ~ hostileMID + polity2 + s_un_reg +  noattackyrs + X_spline1a + X_spline2a + X_spline3a, model="ologit")
# TABLE 4 MODEL 5
z.order5 <- zelig(data=repdata[repdata$polrel==1,], as.ordered(DVordered) ~ hostileMID + polity2 + s_un_reg +  cap_2 + pwrratio + postArt56 + postCW + prgmyrs + NCAnosafetodate + cap_1 + contig + noattackyrs + X_spline1a + X_spline2a + X_spline3a, model="ologit")
# TABLE 4 MODEL 6
z.order6 <- zelig(data=repdata[repdata$polrel==1,], as.ordered(DVordered) ~ hostileMID + postArt56 + postCW + prgmyrs + NCAnosafetodate + noattackyrs + X_spline1a + X_spline2a + X_spline3a, model="ologit")
# TABLE 4 MODEL 7 MY OWN!!
z.order7 <- zelig(data=repdata[repdata$polrel==1,], as.ordered(DVordered) ~ hostileMID + polity2 + s_un_reg +  cap_2 + pwrratio + postArt56 + postCW + prgmyrs + NCAnosafetodate + cap_1 + contig + noconsideryrs + noattackyrs + X_spline1a + X_spline2a + X_spline3a, robust=TRUE,  method="vcovHAC", model="ologit")

summary(z.order1)
summary(z.order2)
summary(z.order3)
summary(z.order4)
summary(z.order5)
summary(z.order6)
summary(z.order7)

?xtable(z.out1,z.out2)
?apsrtable
?glm
# polr.out <- polr(as.ordered(DVordered) ~ hostileMID + polity2 + s_un_reg +  cap_2 + pwrratio + postArt56 + postCW + prgmyrs + NCAnosafetodate + cap_1 + contig + noconsideryrs + noattackyrs + X_spline1a + X_spline2a + X_spline3a,data=repdata,method="logistic")
# beta <- coef(polr.out)
# tau <- polr.out$zeta
# X <- cbind(summary(repdata$hostileMID)[4],summary(repdata$polity2)[4],summary(repdata$s_un_reg)[4],summary(repdata$cap_2)[4],summary(repdata$pwrratio)[4],summary(repdata$postArt56)[4],summary(repdata$postCW)[4],summary(repdata$prgmyrs)[4],summary(repdata$NCAnosafetodate)[4],summary(repdata$cap_1)[4],summary(repdata$contig)[4],summary(repdata$noconsideryrs)[4],summary(repdata$noattackyrs)[4],summary(repdata$X_spline1a)[4],summary(repdata$X_spline2a)[4],summary(repdata$X_spline3a)[4])
# p1 <- plogis(tau[1] - X %*% beta)
# p2 <- plogis(tau[2] - X %*% beta) - plogis(tau[1] - X %*% beta)
# 
# 
# pdf(file = "l4nes1.pdf",width = 6, height = 6, onefile = TRUE,
# family = "Times")
# plot(0:1, p1, type="l", col=1, ylim=c(0,1),las=1,ylab=expression(hat(pi)),xlab="",lwd=2,bty="n",xaxt="n")
# lines(0:6, p2, col=2,lwd=2)
# 
# legend(0, 1, legend=c("P(extremely liberal)", "P(liberal)",
# "P(slightly liberal)", "P(moderate)", "P(slightly conservative)",
# "P(conservative)", "P(extremely conservative)"),
# col=1:7, lty=1,lwd=2)
# Axis(x = 0:6, at = 0:6, side=1,
# labels = c("Strong Dem","WDem","IDem","Indep","IRep","WRep","Strong Rep"))




